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A stochastic version of the Morris-Lecar model is formulated, where 
i i the noise originates from stochastic ion channels. Instead of using a 

\^J Langevin approximation, each population of ion channels [Sodium (Na + ) 

^T' and Potassium (K + )j is modeled as a discrete birth-death process. The 

birth rate is voltage dependent, and the voltage changes 
deterministically between jumps in the number of open ion channels, 
O making the hybrid stochastic process nonlinear. The classic 

I Morris-Lecar model, recovered in the deterministic limit, is an excitable 

system, and transient spikes in voltage are called action potentials. Ion 
channel noise can initiate spontaneous action potentials. In the 
y—( deterministic setting, K + channel dynamics are slow and the fraction of 

^ open K + channels is assumed to be fixed on the time scale of action 

£N| potential initiation. The validity of this assumption in the stochastic 

l/~) setting is examined using a systematic perturbation analysis. We find 

0> that in most physically relevant cases, this assumption breaks down. 
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C^ understand its source 



Any understanding of brain function must include the role of noise. Neural 
networks possess the ability to perform complex computations — taking 
advantage of noise when possible, while still performing reliably. To correctly 
characterize the structure and properties of noise in a network, we must first 



Broadly speaking, a given neuron within a network receives input from two 
main sources of noise. The first is background synaptic activity, such as 
experienced when trying to have a conversation in a crowded restaurant, and 
this noise is a network level phenomenon. The second is spontaneous activity 
due to thermal fluctuations affecting cellular physiology, and this noise is 
intrinsic to each neuron. 
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One source of intrinsic noise is ion channel fluctuations [16]. Sodium (Na + ) 
and Potassium (K + ) ion channels randomly shift between open and closed 
conformations due to the effects of thermal fluctuations, and the rate at which 
channels switch state depends on the membrane voltage. The voltage 
dependent activity of ion channels gives rise to membrane excitability. 
Excitable systems are widely studied in neuroscience. chemical reactions, 
lasers, and climate dynamics [9]. 

Once the voltage crosses a certain threshold, a transient spike in voltage is 
initiated, called an action potential. Ion channel noise can lead to spontaneous 
action potentials (SAPs) , which can have a large effect on network function. If 
SAPs are too frequent, a neuron cannot reliably perform its computational 
role. Hence, ion channel noise imposes a fundamental limit on the density of 
neural tissue. Smaller neurons must function with fewer ion channels, making 
ion channel fluctuations more significant and more likely to cause a SAP. The 
effect of spontaneous activity on the reliability of a neuron can be quantified 
using information theory [14], but the relationship between ion channel noise 
and spontaneous activity remains unresolved. 

In the Morris-Lecar (ML) model of a neuron [10], there is no well-defined 
voltage threshold for initiation of an action potential, but an effective 
threshold can be derived using a fast/slow analysis. In most cases, 
K + channels open and close slowly compared to Na + channels (Ca 2+ channels 
in the original model), and the voltage response to changes in the fraction of 
open Na + channels is so fast that the fraction of open K + channels remains 
relatively constant. However, it is not clear if the fast/slow analysis is valid in 
a stochastic setting. In this letter, we show using a systematic perturbation 
analysis that the fast/slow analysis is not generally valid for initiation of a 
SAP. Furthermore, our analysis of the stochastic system leads to a well 
defined threshold for action potential initiation, allowing for the formulation of 
an exit time problem in an excitable system. 

Deterministic single neuron models, such as the Hodgkin-Huxley model, 
are useful tools for understanding the mechanism of membrane excitability [7]. 
These models assume a large population of ion channels so that their effect on 
membrane conductance can be averaged. As a result, the average fraction of 
open ion channels modulates the effective ion conductance, which in turn 
depends on voltage. The ML model can be understood as a simplified version 
of the Hodgkin-Huxley model and provides a mechanism for the action 
potential. The deterministic ML equation is 
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Koo(v)/n«(v) + wf K (v) + /l ea k(v) + I a pp (1) 
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where v is the membrane voltage, fi{v) = gi{vi — v) determine the ionic 
currents, w is the fraction of open K + channels, and 

Xoo(v) = (1 + tanh(2(7 Na t> + K Na )))/2 is the fraction of open Na + channels, 
assumed to be at quasi-steady state. The steady state for w is 
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Figure 1: (a) Deterministic phase plane dynamics. Thick curves show the null- 
clines: v = as grey and w = as black. Black streamlines represent determin- 
istic trajectories. Green/blue curves represent an action potential trajectory in 
the limit of slow w. (b) Representative stochastic trajectories showing spon- 
taneous action potentials. Stochastic K + and Na + channels: red for voltage 
v(t) and blue for the faction of open K + channels, w(t). Stochastic K + and 
deterministic Na + : orange for voltage v(t) and light blue for w(t). 



Woo(w) = (1 + tanh(2(7 K w + « K )))/2 5 and the time constant 

T w{v) — 2/3 K cosh(7 K u + k k ) is generally assumed to be large so that the w 

dynamics are slow compared to v. We nondimensionalize voltage so that 
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Parameter values are listed in 



v -¥ (v + v cS )/v cS , where 
Appendix A. 

The standard deterministic analysis of an action potential using the ML 
model exploits the separation of time scales. Assume that w is small, and fix 
w = wo constant. The resulting voltage equation is bistable, with a stable 
steady state representing the resting potential, an unstable steady state 
representing the voltage threshold, and a stable steady state representing the 
excited state. In Fig. la, the phase plane dynamics of the system is shown. 
The colored curves that alternate green and blue show the action potential 
trajectory in the limit of slow w. The thick green line at w — wq ~ 0.13 shows 
the initiation phase of the action potential trajectory. Once the excited state is 



reached at the v nullcline (grey curve), the slow w dynamics (blue curve) take 
over. This process continues until the trajectory eventually ends at the fixed 
point. To initiate a SAP, noise must perturb the voltage from the fixed point 
to the starting point of the action potential trajectory. 

Past efforts to understand the relationship between SAP and ion channel 
noise focus on a Langevin (or diffusion) approximation. As a first 
approximation, one can add white noise to a given deterministic equation, 
such as the ML model (1). A more accurate method is to systematically derive 
a Langevin approximation from a more detailed model of ion channel 
fluctuations [3] . However, as recently shown in Ref. [6] , Langevin 
approximations break down when considering metastable dynamics such as 
initiation of a SAP. Moreover, both studies make the same fast/slow 
approximation used to study action potential initiation in the deterministic 
model. 

To consider SAP initiation, one might first add Na channel noise, but 
keep the slow K + channels fixed as in the deterministic analysis. Voltage 
fluctuations caused by stochastic Na + can cause a metastable transition from 
the resting potential over the unstable voltage threshold, initiating a SAP. The 
stochastic initiation of a SAP then reduces to a familiar problem: exit from a 
potential well [6]. Fixing w constant is valid for the deterministic analysis, but 
even if average K + channel activity is slow, how do K + channel fluctuations 
affect SAP initiation? Monte-Carlo simulations of the ML neuron with 
stochastic K + channels (see Fig. lb) show that SAP can be generated by 
K + channel noise alone, without Na + channel noise. If w is not held constant, 
the deterministic system has only one fixed point (at the resting voltage) and 
no longer has a well defined voltage threshold (see Fig. la). Therefore, we are 
immediately faced with a dilemma if we hope to reduce the problem to an exit 
from a potential well. How does one formulate an exit time problem in an 
excitable system with no clearly-defined threshold? If the system is close to a 
saddle-node bifurcation or if w is very slow, one can extrapolate a threshold 
called a ghost separatrix [8]; however, this definition does not hold in general 
and is not an intrinsic property of the excitable stochastic system. 

To answer these questions, we apply recently developed perturbation 
techniques [1,6, 11-13] to the full, discrete stochastic ML model. We find that 
through a stochastic analysis of the ML model, there is a well defined threshold 
for SAP, which does not exist in the deterministic model. Furthermore, we 
find that, in general, w is most likely not constant during initiation of a SAP. 

2 Stochastic Model 

A stochastic version of the ML model is formulated as follows. The voltage 
equation with n = 0, 1, • • • ,N open Na channels and m = 0, 1, • • ■ , M open 

K + channels is 

Ti Til 

V = I(v, m, n) = -/ Na («) + ^/kW + /U(«) + 4pp- (2) 



We assume that each channel is either open or closed and switches between 
each state according to 

O ^± C, i = Na, K, (3) 

where the transition rates are a Na (w) = e 4 (7Na*>+KNa) ; & Na = 1, a K (v) = e ' lKV+KK , 
and b K (v) — e~ 7Kll ~ KK . We assume that the Na + channels open and close 
rapidly, so that l//3 Na <C r m , where r m = C m /g^ is the membrane time 
constant. Taking m and n in (2) to be stochastic birth/death processes, we 
obtain a stochastic hybrid process, which we formulate in terms of the 
probability density function, p(v,m,n,t); it satisfies the differential 
Chapman Kolmogorov (CK) [5] equation, 

d d 

^-p(v, m,n,t) = - — (I(v, m, n)p) + /3 K L K p + /3 Na L Na p, (4) 

at ov 

The jump operators, 

L Na = (E+ - l)fi+ (n\v) + (E r ; - l)n- a (n\v), (5) 

and 

L K = (E+ - m+(m\v) + (E- - 1)0- (m|v), (6) 

govern opening/closing of Na + and K + channels, respectively, with 
E±/(o) = f(a ± 1), S1+ (n|») = n, n^(n|w) = (JV - n)a Na («), 
nj(m|u) = rna K (w), and f2-(m|i>) = (M — m)b K (v). 

The deterministic system (1) is recovered in the limit /3 Na — > oo, M — > oo, 
and we assume that the limit is taken with Am = h^/M fixed. After setting 
w = m/M, the limit yields ^^(f) = o Na (w)/(l + a Na (v)) and 
1000(11) = a K (v) / (b K (v) + a K (u)), which is consistent with the deterministic 
definition. 

The parameter f3 K determines how rapidly the K channels fluctuate. 
Previously, the authors considered slow K + channels [6], with r m /3 K *C 1. 
Here, we assume that v and w change on the same timescale, with 
Tm/3 K = 0(1). The number of open K + channels becomes an external state 
since it shows up as part of the deterministic dynamical system, while the 
number of open Na channels is an internal state since it is averaged out in 
the deterministic limit. 

A perturbation framework has been developed to study metastable activity 
in processes with noisy internal/ external dynamics [1,6,11-13]. Two large 
parameters are present and necessary to reach a deterministic limit, and in 
order to obtain a single small parameter to carry out a systematic 
perturbation expansion, we define e<l such that 

/3- 1 = Xpe, M' 1 = X M e, (7) 



with \fiJT m — 0(1) and A^//r m = 0(1). (We set \p = 1.) Of course, N could 
also be a large parameter, but taking the limit N — > oo, M — > oo yields a 
different deterministic limit than the ML model (1). An alternative to our 
approach would be to take M and N to be large parameters, and this would 
significantly alter the following analysis. However, in the regime where /3 Na , 
M, and N are all large, we expect that either approach yields similar results. 
A more thorough exploration will apear elsewhere. 

We use a WKB perturbation method to derive a uniformly accurate 
approximation of the stationary density [15]. The WKB approximation also 
tells us what path a stochastic trajectory is most likely to follow during a 
metastable transition (i.e., a path of maximum likelihood [4]). First, we 
assume that the stationary solution has the form 



p(v, u>, n) = r(n\v, w) exp 



: <&(v, w) 



(8) 



where <&(v, w) is referred to as the quasipotential and r(n\v, w) is the 
conditional distribution for n given v, w. In the classic problem of exit in a 
double well potential, $ is the double well potential. More broadly, $ is a 
measure of how unlikely it is for a stochastic trajectory to reach a point in 
phase space. After substituting (8) into the right hand side of (4), taking the 
limit t — > oo, and collecting terms in e, we find at leading order, 



L Na +Pv + h(v,w,p w ) 



r(n\v, w) = 0, 



where p v = ^ , and p w — §^ ; and we define 
/3 K 



h(v,w,p u 
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^2(e- jXMP ™ - l)Qi(w\v) 



(9) 



(10) 



j'=± 



where h±(w\v) = 9±(Mw\v)/M. 

In order to solve (9) for $ and r, we first take r to be of the form [1] 



r(n\v, w) 



1 



l + A 



N-r 



A 



l + A 



(11) 



To determine A, we substitute (11) (alternatively, it is simpler to use 

r = A n /n\) into (9) to obtain a consistency expression with two terms: one 

linear in n and one independent of n. From the former we obtain 



A = a Na (u) 



hi 

N 



(p v g(v,w) + h(y,w,p w )), 



(12) 



where g(v,w) — wf K (v) + f\ ea k(v) + I app - After substituting (12) into the 
remaining n-independent term, we obtain the nonlinear scalar PDE for $, 

*^W= ' (13) 



where 

H{v,w,p w ,p v ) 

= -^(1 - x 00 {v))g(v,w)(f Na (v) +g(v,w))pl 
+ (xoo(v)f^(v) +g(v,w))p v 
- -jri 1 - x oo{v))(2g(v, w) + f NB .(v))p v h(v,w,p w ) 

+ h(v,w,p w ) - — (1 -x oo (w))/i(i;,w,p t0 ) 2 , 

which can be solved using the method of characteristics. 

Characteristics are curves (x(t),p(t)) (with x = (v,w) and 
p = (pviPw) = V x $) that satisfy the following Hamiltonian dynamical system, 



x = V p H(x,p), p = -V x H(x,p). (14) 

Note that the deterministic system (1) is recovered by setting p = 0. 
Characteristic projections, x(i), we refer to them as metastable trajectories, 
are paths of maximum likelihood leading away from the deterministic fixed 
point [4]. The action, $(£), satisfying 4>(t) = p(i) • V p "H(x(t), p(t)), is a 
strictly increasing function of t, and the quasipotential is given by 
$>(v,w) = $(£) at the point (v,w) = x(i). Note that along deterministic 
trajectories, p = and $ = 0. We solve (14) using numerical ODE 
integration [11], and we solve (13) on a uniform mesh using an ordered upwind 
finite difference method (OUM) [2]. 



3 Results 

Surrounding the stable fixed point, $ takes the shape of a potential well (see 
Fig. 2), with convex level curves (dashed lines). Once the $ reaches a 
threshold, a point on the level curve looses smoothness forming a caustic. 
Metastable trajectories begin to overlap, and the solution $>(v,w) loses 
uniqueness. Within this region, uniqueness is achieved at each point by 
minimizing the action over all metastable trajectories that pass through that 
point. Using the OUM method, we compute the quasipotential on a uniform 
300 x 300 mesh, for mesh points such that $>(vi,Wj) < $ max ss 0.3. (The 
remaining mesh points are not computed.) The caustic forms an incomplete 
boundary around most of the potential well region. The remaining boundary 
is the curve of constant <&(«, w) = $ c , where $ c w 0.258 is the quasipotential 
at the caustic formation point (see thick dashed line in Fig. 2). We refer to 
this line as the metastable separatrix. 

We identify SAP trajectories as those metastable trajectories that cross the 
separatrix. As shown in Fig. 3, SAP trajectories begin at the fixed point as a 
single trajectory and then fan out just before reaching the metastable 



austic formation point 




Figure 2: Level curves of the quasipotential (dashed curves) and metastable 
trajectories (streamlines). Parameter values are N = M = 10 and Am = 1- 



separatrix (curve marked (S) in Fig. 3). After crossing the separatrix, all of 
the SAP trajectories eventually reach the caustic. Although all SAP 
trajectories are equally likely to reach the separatrix, they are not equally 
likely to reach the caustic. SAP trajectories become less likely the farther from 
the caustic formation point they are when reaching the caustic. Strictly 
speaking, the most probable SAP trajectory strikes the caustic formation 
point, but $ increases by a very small amount in the shaded region of Fig. 3. 
(The maximum difference is |A$| « 0.003, and |A$| /$ c « 0.01.) That is, 
$(t) <C 1 since SAP trajectories are close to deterministic trajectories (black 
streamlines). Hence, the stationary density (8) (after normalization) is nearly 
flat in the shaded region, and all the SAP trajectories in this region are nearly 
equally likely. 

The shaded region represents the most likely experimentally observable 
SAP trajectories; it excludes the small amplitude SAP trajectories that strike 
very close to the caustic formation point and the far less probable SAPs that 
reach the caustic above or behind the potential well region. The SAP 
trajectories that cover the shaded region cross a very small segment of the 
separatrix, the center of which acts much like a saddle point. We refer to this 
point as a metastable saddle. Most of the trajectories that cross the separatrix 
above the saddle are small amplitude SAPs that reach the caustic very close to 
the caustic formation point. Most of the trajectories that cross the separatrix 
bellow the saddle are large amplitude SAPs that are much less likely than 
those in the shaded region. 

Based on the fast/slow analysis of the deterministic system, one might 
expect w to be approximately constant along SAP trajectories if K + channels 




Figure 3: Orange curves are SAP trajectories, shown only up to the point where 
they reach the metastable separatrix (S). The dashed red curve is a SAP that 
crosses S near the metastable saddle (SN); it is not shown after it crosses the 
separatrix because all of the SAP trajectories in the shaded region are visually 
indistinguishable from the dashed red line before they reach S. The shaded region 
contains the most probable SAP trajectories; its lower and upper boundary are 
SAP trajectories that cross S near SN. The SAP trajectories in the shaded region 
are close to deterministic trajectories (black streamlines). Also shown are the 
caustic (C), v nullcline (VN), and w nullcline (WN). Parameters are the same 
as Fig. 2. 



0.5 



0.4 



0.3 



w 



0.2 



0.1 



0.0 




2.0 



Figure 4: SAP trajectories for different values of M and N. Also shown (blue 
curves) are deterministic trajectories. We use Am = l/(eM) with e = 0.1. 



open and close slowly. However, this expectation implies that the SAP 
trajectory is close to a time-reversed deterministic trajectory, which is correct 
only if the stationary density preserves detailed balance. Indeed, the behavior 
of w during a SAP depends strongly on how noisy the K + channels are. Recall 
that K + channel noise is controlled by M since w behaves deterministically in 
the limit M — !• oo. In Fig. 4, several SAP trajectories are shown for different 
values of M and N. Note that the position of the metastable separatrix and 
saddle (not shown) changes for each SAP shown in Fig. 4. If M is very large, 
w behaves deterministically and remains relatively constant during the SAP, 
and in this regime the fast/slow approximation might be valid. On the other 
hand, for most reasonable values of M (i.e., M sa N), the K + conductance 
drops sharply so that v remains below the v — nullcline, and this is also true 
in the limit N — > oo where a SAP is driven entirely by K channel noise 
(dashed curve Fig. 4). 

Fluctuations in the slow recovery dynamics of K channels significantly 
affect the spontaneous activity in the ML model. Since these effects change 
the effective energy barrier leading to a SAP, the spontaneous firing rate will 
depend strongly on the level of K + noise. There are more features of the 
Hodgkin-Huxley model — the sodium conductance has a slow inactivating 
component and ion channels have more complicated multi-state dynamics — to 
be considered in future work. Our analysis uses standard, systematic 
perturbation methods that can be applied to general stochastic excitable 
systems. 

A Parameter values 
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Parameter values are u Na = 120mV, g Na = 4.4mS/cm , vk — — 84mV, 

,9k = 8mS/cm , v\, = — 60mV, g^ = 2mS/cm , C m — 20^iF/cm , 

/3 K = 0.02ms- 1 , I app = 0.06O eff , ip = -0.1, v cS = 52.8mV, 7 », = 1.22/u eff , 

K Na = -1.188 + 1.22w cfF , 7k = 0.8/w Gff , k k = 0.8 + 0.8w c ff. 
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